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Abstract. In this paper, we propose efficient and accurate numerical methods for com¬ 
puting the ground state and dynamics of the dipolar Bose-Einstein condensates util¬ 
ising a newly developed dipole-dipole interaction (DDI) solver that is implemented 
with the non-uniform fast Fourier transform (NUFFT) algorithm. We begin with the 
three-dimensional (3D) Gross-Pitaevskii equation (GPE) with a DDI term and present 
the corresponding two-dimensional (2D) model under a strongly anisotropic confining 
potential. Different from existing methods, the NUFFT based DDI solver removes the 
singularity by adopting the spherical/polar coordinates in Fourier space in 3D/2D, re¬ 
spectively, thus it can achieve spectral accuracy in space and simultaneously maintain 
high efficiency by making full use of FFT and NUFFT whenever it is necessary and / or 
needed. Then, we incorporate this solver into existing successful methods for com¬ 
puting the ground state and dynamics of GPE with a DDI for dipolar BEC. Extensive 
numerical comparisons with existing methods are carried out for computing the DDI, 
ground states and dynamics of the dipolar BEC. Numerical results show that our new 
methods outperform existing methods in terms of both accuracy and efficiency. 
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1 Introduction 


Since its first experimental creation in 1995 (4 20,, 231, the Bose-Einstein condensation 
(BEC) has provided an incredible glimpse into the macroscopic quantum world and 
opened a new era in atomic and molecular physics as well as condensated matter physics. 
It regains vast interests and has been extensively studied both experimentally and theo¬ 
retically |3 17 19 24 34 37 41). At early stage, experiments mainly realize BECs of ultra¬ 
cold atomic gases whose properties are mainly governed by the isotropic and short-range 
interatomic interactions ED- However, recent experimental developments on Feshbach 
resonances [311, on cooling and trapping molecules 38.441 and on precision measure¬ 


ments and control [42.47| allow one to realize BECs of quantum gases with different. 


richer interactions and gain even more interesting properties. In particular, the successful 
realization of BECs of dipolar quantum gases with long-range and anisotropic dipolar in¬ 
teraction, e.g., 52 Cr [261, 164 Dy [35 [ and 168 Er § , has spurred great interests in the unique 
properties of degenerate dipolar quantum gases and stimulated enthusiasm in studying 
both the ground state (7[|8| 29 43 48| and dynamics 113 14 22 27 33 40) of dipolar BECs. 


At temperatures T much smaller than the critical temperature T c , the properties of 
BEC with long-range dipole-dipole interactions (DDI) are well described by the macro¬ 
scopic complex-valued wave function ip = ip(x,t) whose evolution is governed by the 
celebrating three-dimensional (3D) Gross-Pitaevskii equation (GPE) with a DDI term. 
Moreover, the 3D GPE can be reduced to an effective two-dimensional (2D) version if the 
external trapping potential is strongly confined in the z—direction [8,21]. In a unified 
way, the dimensionless GPE with a DDI term in d— dimensions (d = 2 or 3) for modeling 
a dipolar BEC reads as |6][7j[l4,25 48) : 


id t ip(x,t) = 


--V 2 + y(x) + J 6|tp| 2 +AO(x,f) 


ip(x,t), X G lR lf , t> 0, (1.1) 


0(x,f) = (Udip*|i/>| 2 )(x,f), 
ip(x,t = 0) = ip 0 (x) r 


xGK , t > 0, 


xGR rf , 


( 1 . 2 ) 

(1.3) 


where t is time, x= (x,y) r GlR 2 or x= (x,y,z) T E R 3 , * represents the convolution operator 
with respect to spatial variable. The dimensionless constant /3 describes the strength of 
the short-range two-body interactions in a condensate (positive for repulsive interaction, 
and resp. negative for attractive interaction), while V (x) is a given real-valued external 
trapping potential which is determined by the type of system under investigation. In 
most BEC experiments, a harmonic potential is chosen to trap the condensate, i.e.. 


= i f iW+rfy 2 , d= 2, 

2\ 7^ 2 +7jy 2 +7?z 2 , d = 3 


where 7* > 0, 7 y > 0 and 7 Z > 0 are dimensionless constants proportional to the trapping 
frequencies in x-, y- and z-direction, respectively. Moreover, A is a dimensionless constant 
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characterizing the strength of DDI and <I>(x,f) is the long-range dipole interaction whose 
convolution kernel in 3D/2D is given as 


21 281: 


f^dip (*) 


-£(x)-33 r 


47r|x| 


2 (^ n ± n ± n 3 V ±) \2n\x\)' 


77 ffdip (k) — ' 


-Id 


3(n-k) 2 


3[(nj_-k) 2 -ii 2 ||k!| 2 

2||k|| 


d = 3, 
, d = 2, 


(1.5) 


where x,k G and /(k) = J IRrf /(x) e lkx dx is the Fourier transform of /(x). Here, n = 
(hi,H 2 /” 3 ) t is a given unit vector representing the dipole axis, n = (n],ri 2 ) T , 3 n = n-V, 
9nn = 9 n (9n), V_ = (d x ,d lJ ) T , d n± =n ± -V± and d n±n± =d n± (3 n± ). Note that the dipole axis 
can also be different. The dipole kernel L/ c ) j p (x) with two different dipole orientations n 
and m reads as 16,28}|39| 


ffdip (*) — 


— (n-m)^(x) —33 r 


: G1R 3 , 


Jnm V4/r|x| 

-U d ^rn ± -n 3 m 3 \/ 2 L ) ( 2 ^), XG1R 2 , 


( 1 . 6 ) 


where m = ( mi,m 2 ,m 3 ) T is a given unit vector representing the other dipole orientation, 
m — (mi,ni 2 ) T , 9 m± =m^-V_L and =9 n± (9m ± )- We remark here that in most 

physical experiments, the dipoles are polarized at the same direction, i.e., m = n, thus, 
hereafter, we always assume m = n unless specified otherwise. 

The GPE Q-jO) conserves two important quantities: the mass (or normalization) of 
the wave function 

N(t):=\\ip(-,t)\\ 2 := [ \ip(x,t)\ 2 dx = [ \ip(x,0)\ 2 dx = l, t>0, (1.7) 

J R rf J R rf 

and the energy per particle 


E(ip(;t))=[ 

jR d 


^\Vip \ 2 + V(x)\xp \ 2 + ^\xp\* + ^®(x,t)\xp \ 2 


dx=E(xp(- r 0)), t> 0. (1.8) 


The ground state <p g of the GPE ( |1.1[ >-( |1.3[ | is defined as follows: 

</y, = argminE (</>), where S:= {(p(x) \ \\(p\\ 2 := f \cp(x)\ 2 dx = l, E((p) <oo}. (1.9) 

4>eS J Rrf 

Extensive works have been carried out to study the ground state and dynamics of dipolar 
BEC based on the GPE ( |1.1) -(L3). For existing theoretical and numerical studies, we refer 
to 1 7| |17||21||22]|30p2||33 [ and |5,12 I3p8 25]|27 33]|46) , respectively, and references therein. 

To compute the ground state and dynamics of the GPE (|1.1| , one of the key difficulties 
is how to evaluate the nonlocal dipole interaction 0(x,f) ( |1.2| accurately and effectively 
for a given density p=\ip\ 2 . Noticing that 


0(x,f)= / U di Jx-y)p(y r t)dy- 


LT dip (k)p(k,f)e !kx dk, 


/ R rf 


(. 2n) d jRd 


( 1 . 10 ) 
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it is natural to evaluate <5>(x,f) via the standard fast Fourier transform (FFT) using a uni¬ 


form grid on a bounded computational domain 18.39 40.46]. Nevertheless, due to the 
intrinsic singularity/ discontinuity of Udip(k) at the origin k = 0, the so called "numerical 
locking" phenomena occurs, which limits the optimal accuracy on any given computa¬ 
tional domain [ 7|49| . To alleviate this problem, another approach j8)|2l) is to reformulate 
the convolution (1.21 with 3D dipole kernel (|1.5| in terms of the Poisson equation: 


— A«(x,f) = \xp(x,t)\ 2 , lim w(x,t)= 0, xGR 3 , t> 0, 

|x |—>co 


( 1 . 11 ) 


and convolution 0 with 2D dipole kernel ( |1.5| in terms of the fractional Position equa¬ 
tion 

\/—Au(x,t) = | ip (x, t) | 2 , limu(x,f)=0, xGR 2 , t> 0. (1.12) 

|x |—>ca 

Then, the dipolar potential <5>(x,f) can be computed by a differentiation of u(x,t ) as: 


f -|i/>(x,f)| 2 -39 nn »(x,f), xGR 3 , 

<F(x,f) = < t> 0. 

{ -l(9n ± n ± -« 2 V^)w(x,f), XGR 2 , 


(1.13) 


Then in practical computations, the sine pseudospectral method is applied to solve ( 1.11) - 
(1.131 on a truncated rectangular domain Q with homogeneous Dirichlet boundary con¬ 


ditions imposed on 90, and they can be implemented with discrete sine transform (DST) 
efficiently and accurately |8j. By waiving the use of the 0-mode in the Fourier space, 
the sine spectral method significantly improves the accuracy for the dipole interaction 
evaluation. However, due to the polynomial decaying property of u{x,t ) when |x| —l oo, 
a very large computational domain is required in order to achieve satisfactory accuracy. 
This will increase the computational cost and storage significantly for the dipole interac¬ 
tion evaluation and hence for computing the ground state and dynamics of the GPE 0 
Moreover, we shall also remark here that, in most applications, a much smaller domain 
suffices the GPE ( |1.1) simulation because of the exponential decay property of the wave 
function ip(x,t). 

Recently, an accurate and fast algorithm based on the NUFFT algorithm was proposed 
for the evaluation of the dipole interaction in 3D/2D [281. The method also evaluates the 
dipole interaction in the Fourier domain, i.e., via the integral ( 1.10) . Unlike the stan¬ 
dard FFT method, by an adoption of spherical/polar coordinates in the Fourier doma in 
in 3D/2D, the singularity/discontinuity of L7dj p (k) at the origin in the integral (1.101 is 


canceled out by the Jacobian introduced by the coordinates transformation. The integral 
is then discretized by a high-order quadrature and the resulted discrete summation is 
accelerated via the NUFFT algorithm. The algorithm has O(NlogN) complexity with N 
being the total number of unknowns in the physical space and achieves very high accu¬ 
racy for the dipole interaction evaluation. The main objectives of this paper are threefold: 
(i) to compare numerically the newly developed NUFFT based method with the existing 
methods that are based on DST for the evaluation of these nonlocal interactions in terms 
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of the size of the computational domain Q and the mesh size of partitioning Q; (ii) to pro¬ 
pose efficient and accurate numerical methods for the ground state computation and dy¬ 
namics simulation of the GPE with the nonlocal interactions (1.1 l-( |1.2| l by incorporating 
the NUFFT based nonlocal interaction evaluation algorithm into the normalized gradi¬ 
ent flow method and the time-splitting Fourier pseudospectral method, respectively, and 
(iii) to test the performance of the methods and apply them to compute some interesting 
phenomena. 

The paper is organized as follows. In Section 2, we shall briefly review the NUFFT 
based algorithm in |28) for the evaluation of the dipole interaction in 3D/2D. In Section 
3, an efficient and accurate numerical method will be proposed to compute the ground 
state of the GPE ( |1.1| |-( [L2] | by coupling the NUFFT based algorithm for the evaluation of 
the dipole interaction and the discrete normalized gradient flow method. In Section 4, 
we will present an efficient and accurate numerical method for computing the dynamics 
of the GPE 0 -( |1.2[ ) by coupling the NUFFT based algorithm for the evaluation of the 
dipole interaction and the time-splitting Fourier pseudospectral method. Finally, some 
concluding remarks will be drawn in Section 5. 


2 Evaluation of the dipole interaction via NUFFT 

In this section, we will first briefly review the NUFFT based method in [281 for computing 
the dipole interaction in 3D/2D, and then compare this method with the existing DST- 
based method. 


2.1 NUFFT based algorithm 


Due to the external trapping potential, the solution of the GPE (1.11- 0 will decay ex¬ 
ponentially. Thus, without loss of generality, it is reasonable to assume that the den¬ 
sity p(x,t) is smooth and decays rapidly, hence p(k,t) is also smooth and decays fast. 
Therefore, up to any prescribed precision £o (e.g., £o = 10 l2 ), we can respectively choose 
bounded domains V and Br( 0) =: {|k| <R,k£lR‘ f } large enough in the physical space and 
phase space such that the truncation error of p(x,f) and p(k,t) is negligible. Note that the 
convolution only acts on the spatial variable, to simplify our presentation, hereafter we 
omit the temporal variable t and simplify the notation as <f>(x,f)—»<F(x) and p(x,f)—>p(x). 

By truncating the integration domain in ( 1,10} into a Br( 0) and adopting the spheri¬ 
cal/ polar coordinates in 3D/2D in the phase (or Fourier) space, we have 


°( x ) = f e ,kx U dip (k)p(k)dk^ — f e lk ' x U dip (k)p(k)dk 

{2n) a jRd r {2n) a Jb r (q) 

1 f / 0 R /o 27re '' k ' x ! k !^dip(k)p(k)d|k|#, d = 2, 

(2-) | f Q R f 0 n f p TT e' k ' x U d i p (k)\k\ 2 sindp(k)d\k\d6d(p, d = 3. 
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It is easy to see that the singularity/discontinuity of !iai P (k) at the origin is canceled out 
by the Jacobian |k | rf_1 and hence the integrand in the above integral is smooth. High 
order quadratures are then applied to further discretize the above integral and the re¬ 
sulted summation can be efficiently evaluated by the NUFFT [281. The computational 
cost of this algorithm is O(NilogNi) + 0 (N 2 ), where Ni is the total number of equis- 
paced points in the physical space 22 and N 2 is the number of nonequispaced points in 
the phase space Br( 0). Roughly speaking, N 2 is of the same order as N\, however, the 
constant in front of 0 (N 2 ) (e.g., 24 rf for 12-digit accuracy) is much greater than the con¬ 
stant in front of O(NilogNi). This makes the algorithm considerably slower than the 
regular FFT, especially for three dimensional problems. 

To reduce the computational cost, an improved algorithm is also proposed in [281. 
First, by a simple partition of unity, the integral in can be further split into two parts: 


<t ’ (x) “ (2^l { „,^ XQdl '’ (k)?(k) * 

= WfL 0, k " (1 - ^ ( k) > Qdip p C k ),0, e '“ ?<k> dlc 

:= h + I 2 , x6D. (2.2) 


Here, Q = {k=(/c 1 ,...,£: rf ) :r ||k,| < R,l = 1,.. ,,d} is a rectangular domain containing the ball 
Br( 0), the function (j,/(k) is chosen such that it is a C 00 function which decays exponen¬ 
tially fast as |k| —> 00 and the function /y?(k) := (1 — £7 (k ))L/ t jj p ( k ) is smooth for k E R rf . 
With this qd( k), /] can be evaluated via the regular FFT, while b can be computed via the 
NUFFT with a fixed (much fewer) number of irregular points in the phase (or Fourier) 
space (see Figure [lj. Therefore, the interpolation cost in the NUFFT is reduced to 0(1) 
and the overall cost of the algorithm is comparable to that of the regular FFT, with a small 
oversampling factor in front of 0(N\ log/V;). 


2.2 Numerical comparison 


In this subsection, we will show the accuracy and efficiency of the NUFFT based algo¬ 
rithm (referred as NUFFT) for computing dipole interaction and compare them with 
the existing methods that applies ( 1.11| >-( TT3| via the DST (referred as DST). To this end, 
we denote V as the computational domain, 22 ;, as its partition with mesh size h and <5>;, (x) 
as the numerical solution obtained on the domain 22;,. Hereafter, we choose h x = h y = h z 
in 3D and/or h x = h y in 2D and denote them uniformly as h unless stated otherwise. To 
demonstrate the comparison, we define the error function as 


e h :=\\®-® h \ 


I Ol 


(2.3) 


where 


| ;2 IS 


the l 2 - 


norm. 
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Figure 1: Two grids used in the phase (or Fourier) domain in the improved algorithm 
in [281: the regular grid on the left panel is used to compute l \ in ( |2.2| via the regular FFT; 
while the polar grid (confined in a small region centered at the origin) on the right panel 
is used to compute I 2 in ( |2.2| via the NUFFT. Note that the number of points in the polar 
grid is 0(1), thus keeping the interpolation cost in NUFFT minimal. 


Example 2.1. Dipole-dipole interaction in 3D. 

In this example, we take d =3 and choose the source density p(x) =e~ x ~^ j2 with <7>0. 
The 3D dipole interaction with two dipole orientations n and m can be given explicitly 
as 


<f>(x) = — (n-m)p(x) — 39, 


/ 0 2 \J nFxl(r / cr)\ 
V 4 r/a ) 


= — (n-m)p(x) — 3n r G(x)m, 


(2.4) 


where the matrix G(x) = (y ; /(x))? /=1 is given as 


gfli*) 



U^Erf (L 
4 r 3 \cr 


fy+XjXi[- 


r- 

71 



Sfr’yjn 

4r 5 


Erf 



(2-5) 

with Sji the Kronecker delta, x = (xi,X 2 /* 3 ) T and Erf (r) = J Q r e S ~ds the error func¬ 
tion. We choose <7 = 1.4 and compute the potential <t>(x) on a uniform mesh grid, i.e., 
h x = hy = h z on the domain [—L,L] 3 by the DST and NUFFT methods. Table jll shows 
the numerical errors e /, via the DST and NUFFT methods with different dipole axis, 
i.e., n= (0.82778,0.41505,-0.37751 ) T and m = (0.31180,0.93780,-0.15214) T , while Table[| 
presents errors with the same dipole axis, i.e., n = m = (0,0,1 ) r . 

From Tabs. [l]j2j we can clearly observe that: (i) The errors are saturated in the DST 
method as the mesh size h tends smaller and the saturated accuracy decreases linearly 
with respect to the domain size L; (ii) The NUFFT method is spectrally accurate and it 
essentially does not depend on the domain, which implies that a very large bounded 
computational domain is not necessary in practical computations. 
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Table 1: Errors of the 3D dipole interaction with different axis for different h and L. 


NUFFT 

h=2 

h = 1 

h = 1/2 

h = 1/4 

L = 4 

6.004E-01 

6.122E-03 

1.362E-04 

9.823E-05 

L = 8 

6.344E-01 

5.739E-03 

1.189E-09 

6.323E-14 

L = 16 

6.641E-01 

6.054E-03 

1.162E-09 

1.188E-13 

DST 

h = 1 

h = 1/2 

h = 1/4 

h = 1/8 

L = 8 

1.985E-01 

2.022E-01 

2.038E-01 

2.046E-01 

L = 16 

7.135E-02 

7.172E-02 

7.200E-02 

7.214E-02 

L = 32 

2.622E-02 

2.544E-02 

2.549E-02 

2.552E-02 


Table 2: Errors for the 3D dipole interaction with the same axis for different h and L. 


NUFFT 

h = 2 

h = 1 

ft = l/2 

/z = 1/4 

L = 4 

1.118E-01 

3.454E-04 

1.335E-04 

1.029E-04 

L = 8 

5.281E-02 

3.428E-04 

9.834E-12 

1.601E-14 

L = 16 

5.202E-02 

3.551E-04 

1.143E-11 

8.089E-15 

DST 

h = 1 

h = 1/2 

h = 1/4 

h = 1/8 

L = 8 

6.919E-02 

7.720E-02 

8.124E-02 

8.327E-02 

L = 16 

2.709E-02 

2.853E-02 

2.925E-02 

2.961E-02 

L = 32 

1.008E-02 

1.033E-02 

1.046E-02 

1.052E-02 


Example 2.2. Dipole-dipole interaction in 2D. 


Here, we take d = 2 and choose the source density as p(x) = e with n>0. The 

exact 2D dipole interaction with two dipole orientations n and mi can be given as 


O(x) = 


3y/ Tie 
4 a 

1+2 r 
2s~ 


( n J 


h (r 


™±)(k(r)-h(r))~ 
?>^/nn 3 m 3 re 


2(x-ni)(x-m j) 


k(r) 


h(r)-h(r)- 


k{r) 

2 r J' 


( 2 . 6 ) 


where r = ^ T 


, Iq and / ] are the modified Bessel functions of order 0 and 1, respectively HI- 
Here, we choose <j=\J 1.3 and dipole axis as = (0,—0.896) r and m^ = (0,—0.52476) T 
(corresponding to n = (0,—0.896,0.44404 ) T and m= (0,—0.52476,0.85125 ) T in 3D). Table 
[3] shows the errors t’/, via the DST and NUFFT methods for different domain size L and 
mesh size h. 

From Tab. |3j we can clearly observe that: (i) The errors are saturated in the DST 
method as the mesh size h tends smaller and the saturated accuracy decreases linearly 
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with respect to the domain size L; (ii) The NUFFT method is spectrally accurate and it 
essentially does not depend on the domain if it is adequately large. 

Table 3: Errors of the 2D dipole interaction by different methods with h on [— L,L] 2 . 


NUFFT 

h=2 

h = 1 

h = 1/2 

h = 1/4 

L = 4 

11.96 

6.444E-01 

5.251E-06 

7.343E-06 

L = 8 

1.279 

1.611E-02 

4.039E-07 

4.720E-14 

L = 16 

3.289E-01 

1.631E-02 

4.226E-08 

3.489E-14 

DST 

h = 2 

h = 1 

h = 1/2 

h = 1/4 

L = 8 

3.200E-01 

1.944E-02 

1.145E-02 

1.208E-02 

L = 16 

3.264E-01 

1.660E-02 

2.971E-03 

3.048E-03 

L = 32 

3.281E-01 

1.636E-02 

7.590E-04 

7.686E-04 


Example 2.3. Dipole-dipole interaction for anisotropic densities. 

In this example, we consider the dipole-dipole interaction for anisotropic densities 
which are localized in one or two spatial directions. As stated in the introduction, the 
3D/2D dipole interaction ( 1.13) can be solved analytically via the Poisson equation ( 1.11) 
and the fractional Poisson equation ( |1.12 l in 3D and 2D, respectively. Therefore, the 
dipole-dipole interaction can be obtained analytically via the solution of the Poisson/ 
fractional Poisson potential, followed by differentiation. 

i * 2 y 2 

The 2D case. For an anisotropic density p(x,y ) = 4 with a small parameter 

0 < £ < 1, the 2D Coulomb potential ( 1.12) is given analytically [16] as: 


u(x,y) = 


2n\fn Jo 


G(x,y,s)ds, G(x,y,s) = 


exp( — 


4(l+s 2 ) 4 (s 2 +£ 2 ) 


) 


\/s 2 +1 \/s 2 + £ 


Then, the 2D DDI can be obtained by differentiating the integrand G in 
convenience of readers, it can be evaluated explicitly as: 


(2.7) 


For the 


3 f°° 

0 (x,t/) = -^ 3 7 2 J o ((nim-L -n 3 m 3 )G xx +(n 2 m 2 -n 3 m 3 )G !/ y+(n 1 m 2 +n 2 m 1 )G xy )ds. (2.8) 


Similarly as [161, to numerically evaluate ( |2.8) , we first split it into two integrals and 
reformulate the one with infinite interval into an equivalent integral with finite interval 
by a change of variable. Then, we apply high order Gauss-Kronrod quadrature to each 
integral to get reference solutions. Here we omit details for brevity. With this way, we 
can obtain the 'exact' 2D DDI with the given density p(x,y). 

As discussed in [161, the 2D Coulomb interaction can be well-resolved by the NUFFT 
method on a heterogenous rectangle T> e =[—L,L]x [— eL,eL\. The DST method, best suited 
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for solving the PDEs with homogeneous Dirichlet boundary condition on a rectangular 
domain, fails to produce even a satisfactory result on T> e , mainly because the potential 
does not decay fast enough. Actually, the DST method can give reasonably accurate re¬ 
sults on a square V= [—L,L] 2 due to that the homogeneous Dirichlet boundary condition 
doesn't bring significant error, however, one needs to resolve the anisotropic density with 
h y = eh x . Then the computational and storage costs for the DST method will correspond¬ 
ingly scale linearly as a function of 1 / e. Here we adapt the similar strategy for the choice 
of computational domains for the NUFFT and DST methods to compute the DDI. Table 
[4] presents errors of the 2D dipole interaction for anisotropic densities by NUFFT on 'D c 
and DST on V with h x = 1/8 ,h y = eh x for the same n,m as in the Example 

Table 4: Errors of the 2D dipole interaction for anisotropic densities solved on V £ and V 
for the NUFFT and DST methods, respectively, with h x = 1/8. 


2.2 


NUFFT 

£ = 1 

£ = 1/2 

£ = 1/4 

£ = 1/8 

£ = 1/16 

L = 8 

3.456E-07 

4.167E-07 

3.984E-07 

3.207E-07 

2.864E-07 

L = 16 

1.005E-12 

7.531E-15 

5.119E-15 

4.108E-15 

3.720E-15 

L = 32 

2.241E-12 

6.856E-15 

4.913E-15 

4.072E-15 

3.855E-15 

DST 

£ = 1 

£ = 1/2 

£ = 1/4 

£ = 1/8 

£ = 1/16 

L = 8 

4.014E-02 

3.182E-02 

1.711E-02 

6.276E-03 

2.181E-03 

L = 16 

9.604E-03 

7.534E-03 

4.035E-03 

1.479E-03 

5.137E-04 

L = 32 

2.386E-03 

1.868E-03 

9.995E-04 

3.661E-04 

1.272E-04 


As for the 3D case, there are two typical kinds of anisotropic densities, that is, densi¬ 
ties that are strongly localized in one or two directions. The first typical kind of anisotropic 
density is localized in one direction. For example, choose the density as p{x,y,z) = 


87 ly/m 


x 2 +y 2 _ 

e and its corresponding 3D Coulomb potential (|1 .111 is given as: 


u{x,y,z ) : 


87T 3/2 . 


£ 4(l+s) £ 4(s+e 2 ) 


(1 + s) \/s + £' 


-As. 


(2.9) 


The second kind of anisotropic density is localized in two directions. For example, the 

-1 _ * 2 +!/ 2 z 2 

density is taken as p(x,y,z) = 4,2 and the corresponding 3D Coulomb po¬ 

tential (|1.11| is given analytically as: 


u(x,y,z) = „ / e 4(s+ f 2 ) e 415 +t) -— 

8tt 3/2 Jo yT+s (s+e 2 ) 


ds. 


( 2 . 10 ) 


Similarly as in the 2D case, the DDI in 3D can be obtained by differentiating the integrand 
in \2.9) and ( 2.1 0[ ). They can be evaluated numerically in a similar way, which can be 
viewed as the 'exact' solution. For brevity, we omit the formulas and relevant details. 
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To numerically compute the 3D dipole interaction by the NUFFT and DST methods, 
we shall adopt the meshing strategy, i.e., h y = h x and h z = eh x for densities localized only 
in z-direction and h y = h x = eh z for densities localized in x,y directions. Similarly, the 
NUFFT method is applied on a heterogenous cube T> e = [—L,L] 2 x [— eL,eL] or [—eL,eL] 2 x 
[— L,L\. The DST method is used on V = -L,L] 3 so that the homogeneous Dirichlet 

boundary condition doesn't bring significant error. Correspondingly, the computational 
and storage costs of the DST method will scale linearly as a function of 1 /e (the first kind) 
or 1/e 2 (the second kind). 

To show the accuracy performance of both methods, we take the first kind density as 
the test function Here we take the same dipole axis n = m = (0,0,1 ) T for simplicity. 
Table [5] presents errors of the 3D dipole interaction for anisotropic densities localized in 
the z-direction by NUFFT on V £ and DST on V with h x — 1/4 and h z = eh x . 


Table 5: Errors of the 3D dipole interaction ( |2.9| by NUFFT on V £ and DST on V with 
h x = 1/4 and n = m = (0,0,l) r . 


NUFFT 

£ — 1 

£ = 1/2 

£ = 1/4 

£ = 1/8 

£ = 1/16 

L = 8 

3.004E-08 

2.581E-08 

2.307E-08 

1.988E-08 

1.578E-08 

L = 16 

1.598E-14 

7.590E-15 

4.590E-15 

2.184E-15 

1.193E-15 

DST 

£ = 1 

£ = 1/2 

£ = 1/4 

£ = 1/8 

£ = 1/16 

L = 8 

1.409E-01 

7.667E-02 

4.308E-02 

2.633E-02 

1.716E-02 

L = 16 

5.003E-02 

2.754E-02 

1.548E-02 

9.453E-03 

6.159E-03 

L = 32 

1.786E-02 

9.836E-03 

5.522E-03 

3.370E-03 

2.195E-03 


From Tabs. |4]j5j we can conclude: (1) the NUFFT can evaluate accurately the 2D and 
3D dipole interaction with anisotropic densities. (2) The DST method can still capture 
satisfactory results if the computational domain is large enough, however, the computa¬ 
tional and storage costs increase when the heterogeneity of the density increases, which 
makes it less applicable, especially in 3D simulation. 


3 Ground state computation 

In this section, we propose an efficient and accurate numerical method for computing 
the ground state by combining the normalized gradient flow which is discretized by the 
backward Euler Fourier pseudospectral method and the NUFFT nonlocal DDI interaction 
solver. We shall refer to this new method as GF-NUFFT hereafter. The spatial spectral 
accuracy is investigated in details, the virial identity is verified numerically, with com¬ 
parison to some existing results in [8|, to show the advantage of the GF-NUFFT method 
in term of accuracy. 
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3.1 A numerical method via the NUFFT 

Let At > 0 be the time step and denote t n = riAt for n = 0,1,2,... . Many efficient and 
accurate numerical methods have been proposed for computing the ground state of the 
GPE ||8]- fT0||49| . One of the most simple and successful method is by solving the following 
gradient flow with discretized normalization (GFDN): 


dt<p{x,t) = 


-V 2 -V(x)-p\xp\ 7 


0(x,f) = (U dip *|</>| 2 )(x,f), 

(piX'tn+l) 


<p{*,tn+ 1 ) : 


iw*,f„ + i)ir 


AO(x,f) 


<p(x,t), xGR rf , 


x G R‘ f , 


xGr, 


tn^t < fji+l, 
— tji+lr 

n> 0, 


(3.1) 

(3.2) 

(3.3) 


with the initial data 


0(x,O) = 0 o (x), xGE'', 


with ||^o|| 2 := [ \(po{x)\ 2 dx = l. 


(3.4) 


Let 0"(x) and <3>' ! (x) be the numerical approximations of < p(x,t n ) and <E>(x,f„), respec¬ 
tively, for n > 0. The above GFDN is usually discretized in time via the backward Euler 


method [|8}[49[ 


(x) — <p n (x) 

At 

0"(x) = (U dip *|</>' ! | 2 )(x), 

Hr+i( x )= ^ (1) W 

* (x) n^)(x)ir 


^V 2 -y(x)- J 6|^"| 2 -AO"(x) 


tf> (1) (x), xGR d , 


xGR d , 


x G R , n > 0. 


(3.5) 

(3.6) 

(3.7) 


As it is known, the ground state decays exponentially fast due to the trapping poten¬ 
tial, therefore, in practical computations, we shall first truncate the whole space to a 
bounded domain V and impose periodic boundary conditions. It is worthwhile to point 
out that the dipole interaction is originally defined by convolution, therefore it does not 
require any boundary condition. Then, the equation ( |3.5| > is discretized in space via 
the Fourier pseudospectral method and the dipole interaction 0"(x) in ( ]3.6| > is evalu¬ 
ated by the NUFFT solver. The initial guess (po{x ) is usually chosen as a positive func¬ 
tion, e.g., a Gaussian, and the ground state (pg(x) is obtained numerically as ([>” (x) once 

^ — ' X ^'” G £o is satisfied, where £o is the desired accuracy. The details are omitted 

here for brevity. 
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3.2 Numerical results 

In order to study the spatial accuracy of the GF-NUFFT method for computing the ground 
state, we denote <h ? (x) = (L/j ip * \(p g \ 2 )(x) and introduce the error functions 

* _ II^W-^WIIp h _ ll°^( x )-°g( x )ll / 2 
IlfeMllp ' ll*,(x)|| P ' 

where (p g and <E>g are obtained numerically by a numerical method with the mesh size h. 
Additionally, we split the energy functional into four parts 


E((p) — Ekin(^) + Epot (</>) +hint (<^) + £(%(</>)> 

where the kinetic energy E k j n (<^), the potential energy E pot ((p), the interaction energy 
Emt{(p), and the dipole interaction energy E L y ip (ip) are defined as 


£kin (<p) = 2 Jr* I V(p(x)j 2 dx, E pot (<p) = f Rd V (x) | (p(x) | 2 dx, 
hint (cp) = § f Rd l<p(x) I 4 dx, Edip(^) = 2 f Rd ®(x)\(p(x) | 2 dx. 


respectively Moreover, the chemical potential can be reformulated as }i((p) = E(cp) + 
Emt((p) +Edip(<^). Furthermore, if the external potential V(x) in (1.1) is taken as the har¬ 
monic potential ( |1.4[ ) 015 [16 36 [, the energies of the ground state satisfy the following 
virial identity 

0 = I ■ = 2Ekm( ( | , g) —2Ep 0 t((pg) + 3Ei nt (<|>g) + 3Edi p (^g). (3-8) 


We denote I h as an approximation of I when <p g and Og are replaced by cp 1 ', and in (3.8 . 

In our computations, the ground state (p g is reached numerically when —— 
eo with £o = 10“ 10 . The initial data (po(x) is chosen as a Gaussian and the time step is taken 
as Af = 10 2 . In the comparisons, the "exact" solution <pg(x) was obtained numerically 
via the GF-NUFFT method on a large enough domain Q with small enough mesh size h 
and the same time step A t = 10 2 . 

Accuracy confirmation. To show the accuracy of the GF-NUFFT, we take d = 3, 
n= (0,0,1) T , p = 200 and V(x) = 2(x 2 +y 2 +z 2 /4). Table |6] presents errors of the ground 
states and the corresponding dipole interactions computed on a fixed domain [—8,8] 3 
with different mesh sizes and A. From this Table, we can observe clearly the spectral 
convergence in space of the GF-NUFFT method. 

Virial identity. Here we take the same physical parameters as used in [8j (cf. Table 
3), i.e., d = 3, /3 = 207.16 and V(x) = ^{x 2 +y 2 +z 2 /4). We compute the ground state on 
a larger domain, i.e., [—12,12] 3 , with a coarser mesh size h x = h y = h z = 1/4. Different 
energies of the ground state and related quantities are shown in Table [7j Compared with 
Table 3 in [8| where the identity is only accurate up to 3 significant digits, our results by 
the GF-NUFFT method agree quite well with the identity, up to 9 significant digits. 
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Table 6: Errors of the ground states and the dipole interaction obtained by the GF-NUFFT 
method for n= (0,0,1 ) T and /l = 200 with different mesh sizes h and A. 


GF-NUFFT 

h = 2 

h = 1 

h = 1/2 

7i = 1 / 4 

h = 1/8 


A = 

-100 

1.783E-02 

3.102E-03 

3.463E-05 

3.652E-09 

4.133E-12 

ph 

A = 

100 

1.263E-02 

2.717E-03 

3.599E-05 

6.183E-09 

2.841E-12 

• o 

A = 

200 

1.670E-02 

3.049E-03 

8.897E-05 

5.364E-08 

3.871E-12 


A = 

-100 

2.810E-02 

3.683E-03 

1.842E-05 

1.555E-09 

8.132E-12 

ph 

A = 

100 

2.385E-02 

4.932E-03 

9.445E-05 

1.996E-08 

2.750E-12 

8 

A = 

200 

1.406E-02 

5.681E-03 

2.424E-04 

1.921E-07 

3.121E-12 


Table 7: Different energies of the ground state and I 1 ' for the 3D dipolar BEC with /3 = 
207.16 for different A. 


A 


l l s 

F 8 

E kin 

E 8 

^pot 

E 8 

^ int 

E 8 

c di P 

I h 

-103.58 

2.9584 

3.9301 

0.26466 

1.7221 

0.83892 

0.13273 

6.6214E-10 

-51.79 

2.8841 

3.8187 

0.27379 

1.6757 

0.85255 

0.082056 

5.7861E-10 

0 

2.7943 

3.6830 

0.28621 

1.6193 

0.88875 

0.0000 

5.0929E-10 

51.79 

2.6875 

3.5201 

0.30303 

1.5519 

0.94903 

-0.11646 

4.5134E-10 

103.58 

2.5593 

3.3213 

0.32704 

1.4701 

1.0451 

-0.28304 

3.6672E-10 

155.37 

2.3998 

3.0674 

0.36538 

1.3668 

1.2105 

-0.54290 

2.3288E-10 

207.16 

2.1838 

2.7011 

0.44525 

1.2212 

1.5749 

-1.0576 

-1.6697E-10 


4 Dynamics simulation 

In this section, instead of solving the GPE we consider a more general GPE in 

d-dimensions (d = 2,3) with both the damping term and time (in-)dependent DDI: 


id t xp(\,t) 


-yV 2 + y(x) + ^|i/;| 2c7 +AO(x,f)-//(|ip| 




0(x,f) = (Lr dip *|j/>| 2 )(x,t), xGR rf , f>0, 
ip(x,t=0) = xp 0 (x). 


(4.1) 

(4.2) 

(4.3) 


Here, cr > 0 corresponds to the type of the nonlinearity (a = 1 represents to the cubic 
nonlinearity, and resp., a = 2 to a quintic nonlinearity). f(p) > 0 for p = \ if)\ 2 > 0 is a real¬ 
valued monotonically increasing function that represents the type of damping. In BEC, 
when f(p) =0, ( |4.1| reduces to the usual GPE (JtTJ without damping effect, while a linear 
damping term f(p ) = 3 with 5 > 0 represents inelastic collisions with the background gas. 
In addition, the cubic damping f(p) = 3\p with Aj > 0 describes two-body loss, a quintic 
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damping term of the form f(p) = <$ 2 p 2 with 82 > 0 corresponds to the three-body loss, and 
their combination f(p) = p + d 2 p 2 takes both the two and three-body loss into account. 
Furthermore, the kernel of the dipole interaction may be time (in-) dependent, which is 
defined as 


f^dip 


3 m(f) -n(f) — 3(x-m(f))(x-n(f))/|x | 2 
47T |x| 3 

— ( m (f)-n(f))(5(x) — 3d m(t)n(t) ( 4 ^)' xGK 3 , 


(4.4) 


where m(f) = (mi(f) r m. 2 {t),m. 2 ,{t)) T and n(f) = (ni(f),W 2 (f), ri 3 (t)) T G 1R 3 are two given 
time (in-)dependent unit vectors, representing the two dipole orientations. The energy is 
modified as: 


m 



|Vi/?| 2 + V(x)|i/;| 2 G 


s+i 


\ip\ 2 ( s+1 '> 


^0(x,f)|ip| 2 


A r t / 

— 7T / ( ds Fldip * I 


2 Jo 


\xp(x,s)\ 2 ds 


dx. 


which satisfies the following dynamical law: 


2 . 


E‘ f 


f(\xp\ 2 )lm(ipd t ip)d 


x. 


(4.5) 


(4.6) 


where ip denotes the complex conjugate of ip. The total mass N(t) (1.7' is dissipated as: 

d 


-N(t)=-ij^m\ 2 M 2 dx. 


(4.7) 


We will present an accurate and efficient numerical method for simulating the dy¬ 
namics of the GPE g-g. The method incorporates the NUFFT solver for the eval¬ 
uation of the nonlocal dipole interaction and the time-splitting Fourier pseudospectral 
discretization for the GPE ( |4.1| . 


4.1 Numerical method 

In practical computation, we first truncate the problem ( |4.1| |-( |4.3[ ) into a bounded com¬ 
putational domain T>= [L X ,R X ] x [L y ,Ry] x [L Z ,R Z ] if d = 3, or V= [L X ,R X ] x [L y ,R y ] if d = 2. 
From t = t n to t = t n+ the GPE ( |4.1j ) will be solved in two steps. One first solves 

id t ip(x,t) = --V 2 xp(x r t), xGP, t n <t<t n+ 1 , (4.8) 

with periodic boundary condition on the boundary dV for a time step of length At, fol¬ 
lowed by solving 

id t ip(x,t) = [V(x) + J 6|!/>| 2l7 + AO(x,f)-z/(|i/?| 2 )] ip(x,t), xEV, t n < t < t n+1/ (4.9) 

0(x,t) = (U dip *\xp\ 2 )(x,t), xeV, f„<f<f„+ 1 , (4.10) 
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for the same time step. The linear subproblem ( |4.8| > will be discretized in space by the 
Fourier pseudospectral method and integrated in time exactly in the phase space, while 
the nonlinear subproblem (|4.9|>-(|4.10]> can be integrated exactly, one can refer to [8.11.13 


141 for details. To simplify the presentation, we will only present the scheme for the 3D 
case. As for the 2D case, one can modify the algorithm straightforward. 

Let L, M, N be even positive integers, choose h x = Rx y Lx , h y = Ry M L| ' and h z = Rz N L: as 
the spatial mesh sizes in x-, y-, and z- directions, respectively. Define the index and grid 
points sets as 


Tlmn = {(/,fc,m)|0</<L, 0 <k<M, 0<m<N}, 
m )< m i N ^ ^ N it 

Gxyz = { {xi,yk>Zni) = - {F x -\-jh x ,L y -\-khy,L, z -\-mh z ), (l,k,m) eTlmn} ■ 

Define the functions 

(p, q ,r) e f LMN , 

with 


v-y 


2 pn y 

Tyr Lv 


1 TT_ u z. 


2rn 


R y Ly 


Ry — L 


-, ( p,q,r)eT L MN ■ 


Let ipf km be the approximation of ^(x/,y fc ,z m ,f n ) for ( l,k,m ) E Tlmn arid n >0 and denote 
ip n be the solution vector at time t = t n with components ( l,k,m ) eTlmn}- Taking 

the initial data as xp^ km = ipo(xi,yi c ,z m ) for (. l,k,m ) E Tlmn, for n >0, a second-order time 
splitting Fourier pseudospectral (TSFP) method to solve the GPE ( |4.1[ )-( |4.3| reads as: 

, L/2-1 M/2-1 N/2-1 , A , r „ „„ _ 

fill = E E E 

p=—L/2q=—M/2r=—N/2 

W2« = Wto« ex p{- I, [ A * vr ( x )+ H (IV’/2«l 2 ' A 0+G(|^ (1) | 2 ,t n ,t" +1 )(ac/,yfc,z m )]} 

xexp{—F(|i/i|^J 2 ,Af)}, (l,Km) ETlmn, 

L/2-1 M/2-1 N/2-1 ,., r „ „„ 

ffi 1 = E E E e - , m )2 +« )2 +« ,2 ](^>) wr w*ft*,,y b zj. (4.11) 

p=—L/2q=—M/2r=—N/2 

Here, (tp n ) pqr and (^ (2) ) p(?r are the discrete Fourier transform coefficients of the vectors 
xp n and (/r 2i , respectively, and the functions H(cp,s), G(cp,s,s\) and F(<p,s) are defined as: 


H{<p,s) = P[h(cp,r)] a dT, F(f,s) = J f{h(q>,r))dT, 

fS i 

G(^,s,si)(x)=A J (H dip (-,T)*/z(<p(-,T-s)) (x)dr, 


(4.12) 


(4.13) 
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h(f,s) = 


*00 = /; 


(4.14) 


with 

g~\g(<p)~ 2s), <p>0, s >0, 

0, <p = 0, s>0, J sf(s) 

For a given damping function f(s), in general, g 1 (s) and thus h{cp,s) may not have 
explicit expressions. In practical computation, one could solve h(<p,s) numerically from 
an auxiliary ODE, and then evaluate ( |4.12 )-(4.131 via a numerical quadrature. For details, 
one can refer to Remark 2.1 in ('111. However, if the dipole axis is time independent, i.e., 
lidip(x,f) = Lfdi p (x,f = 0) =: H° ip (x), for those damping terms that are frequently used in 
the physics literatures, the functions H, F and G can be integrated analytically. For the 
convince of the reader, we list them here briefly (lT( : 

• Case I. f(p) = 0, i.e., no damping term, we have 

H((p,s) = pf <r s, F(<p,s)= 0, G(<p,s,si)(x)=A(si-s) (x). 


• Case II. f(p)=5> 0, i.e., the linear damping, we have 

H(cp,s) = F^(l-e~ 2Scrs '), F(tp,s)= 5s, 

G(f,s,»i)(x) = 4 (uSi P *f>) (*)• 


• Case III. f(p) = Sp q with S,q > 0, which corresponds to two (q = 1) or three (q = 2) 
body loss of particles, we have 

F(cp,s ) = -*-ln(l+2(7&<p' 7 ). 


H(cp,s) = 


2^1n(l+2 qSsqfl), 


pep 1 ' 1 
2S(q—a) 


(l+2q5s(p CJ ) 1 ~ a/ci — 1 


if cr = q, 
if u^q, 


G(cp,s,s 1 )(x) = \U° dip 



3jln(l+2£(si-s)<p), 

(l+2 9 ,5(s 1 -s)<p?) 1 ^ 1/? -l 
2S(q—l)cp^ 1 ' 


if £7 = 1, 
if <7^1. 


The function G is evaluated by the algorithm via the NUFFT as discussed in previous 
sections, and this method for discretizing the GPE (|4.1[l-||4~3|) is referred as TS-NUFFT. 


4.2 Test of the accuracy 


In this section, we first test the accuracy of our numerical method for computing the 
dynamics of the dipolar BEC. To demonstrate the results, we define the following error 
function 


h,M 

e <P 


_ ||t/(x,f n )-^ Af (x)|| F 

IIiKx,0IIp 


n> 0, 


(4.15) 
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where || • || p represents the l 2 norm, Af (x) is the numerical approximation of tp(x,t = t n ) 
obtained by the TS-NUFFT method (4.111 with mesh size h and time step At. In this 
subsection, all examples are carried out for dipolar BEC without damping effect, i.e., 
f(p)= 0 in the GPE (4.11. Moreover, the computational domain 'D, the trapping potential 
V (x) and the initial data i/’o(x) are respectively chosen as 


V=l- 2 ^ d ,^- 4 ] d , V(x) = !i, lfo(x) = 4 ? 


. M 


e 2 , xeV with d = 2 or 3. (4.16) 


Furthermore, the dipole orientations are chosen as n = m = (0,0,1 ) T in 3D and n^ = m = 
(1,0) T in 2D, respectively. 


Spatial errors (upper parts) e^' Af °(f) and temporal errors (lower parts) e 


Table 8 

t = 0.28 for the dynamics of the 3D GPE (4.1 1 with different /3 and A = j ■ 


h 0 ,At 


(0 at 



h = 1/2 

II 

5 I 

ft = 1/8 

ft = 1/16 

(6 = 2 

3.999E-03 

1.612E-05 

1.601E-11 

3.049E-12 

j8 = 10 

1.773E-02 

2.581E-04 

8.899E-09 

3.133E-12 

j8 = 50 

8.074E-02 

8.186E-03 

2.460E-05 

2.304E-11 

4°' A£ (0 

> 

if 

o 

o 

o 

00 

At/2 

Af/4 

At/8 

/3 = 2 

2.983E-06 

7.454E-07 

1.860E-07 

4.615E-08 

rate 

- 

2.001 

2.003 

2.011 

/3 = 10 

8.151E-06 

2.036E-06 

5.081E-07 

1.261E-07 

rate 

- 

2.001 

2.003 

2.011 

j8 = 50 

8.427E-05 

2.105E-05 

5.251E-06 

1.303E-06 

rate 

- 

2.001 

2.003 

2.011 


Example 4.1. Numerical accuracy verification in 3D. 

Here d=3 and the "exact" solution ip(x,t) is obtained numerically via the TS-NUFFT 
method on domain T> with very small mesh size h=ho: = j^ and time step Af=Afo:=10~ 4 . 
Tablejijlists the spatial discretization errors e^' Af ° ( t ) and the temporal discretization errors 
ejjj ),A, (T) as well as the convergence rate at time t = 0.28 with different mesh size h and 
different time step At, for different /3 and A= j. 

Example 4.2. Numerical accuracy verification in 2D. 

Here d=2 and the "exact" solution ip(x,t) is obtained numerically via the TS-NUFFT 
method on domain T> with very small mesh size h=hy. = ^ and time step Af=Afo:=10~ 4 . 
Table 9 shows the spatial discretization errors e^ /Af °(t) and the temporal discretization 
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errors e^°' At (t) as well as the convergence rate at time t = 1.0 with different mesh size h 
and different time step At, for different /I and A = J). 


Table 9: Spatial errors (upper parts) e^ At ° ( t ) and temporal errors (lower parts) e’p 


llQ ,Af 


t = 1.0 for the dynamics of the 2D GPE (4.1' with different (1 and A = 


(0 at 


ei At0 (t) 

h = 1/2 

h = 1/4 

OO 

T-1 

II 

h = 1/16 

P=2 

5.715E-05 

6.193E-11 

1.120E-11 

1.124E-11 

16 = 10 

1.894E-03 

6.616E-08 

1.354E-11 

1.679E-11 

j6 = 50 

7.265E-02 

2.987E-04 

4.987E-10 

2.852E-11 


Af = 0.01 

At/2 

Af/4 

At/8 

P = 2 

9.011E-06 

2.252E-06 

5.623E-07 

1.399E-0 7 

rate 

- 

2.001 

2.002 

2.007 

II 

h-* 

O 

2.293E-05 

5.728E-06 

1.430E-06 

3.558E-07 

rate 

- 

2.001 

2.002 

2.007 

j6 = 50 

2.453E-04 

6.122E-05 

1.528E-05 

3.802E-06 

rate 

- 

2.003 

2.002 

2.007 


From Tabs. [8]j9j we can see that the TS-NUFFT method ( 4.11| is spectrally accurate in 
space and second order accurate in time for computing the dynamics of dipolar BEC. 


4.3 Applications 


In this section, we apply the TS-NUFFT method ( |4.11 1 to study some interesting phe¬ 
nomena, such as the dynamics of a BEC with time-dependent dipole orientations and the 
collapse and explosion of a dipolar BEC with attractive interaction and damping terms. 

Example 4.3. Dynamics of a BEC with rotating dipole orientations. 


Here d=3 and we consider the GPE (|4.1 } without damping term, i.e., f(p) = 0. The 

| 1 2 *-* 

trapping potential is chosen as V(x) = ^- and the i nitia l data in (4.31 is chosen as ipo (x) = 
0 gs (x), where <^ gs (x) is the ground state of the GPE (4.1 1 with/(p)=0 and n=m=(0,0,l) T , 


j6 = 103.58 and A = 82.864, which is computed numerically via the numerical method 
presented in the previous section. The computational domain and mesh size are chosen 
asV= [—8,8] 3 and h x = h y = h z = |, respectively. Then we tune the dipole orientations as 


, \ , t t 

n(t) = sm-,0,cos- 
' 5 5 


f>0. 


(4.17) 


and study the dynamics of the BEC in two cases: 
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Case 1. tune the dipole orientation as in ( 4.17| and keep all the other parameters 
unchanged. 

Case 2. tune the dipole orientation as in ( 4.17| , perturb the trapping potential by 
setting 7 * =2 and keep all the other parameters unchanged. 


Figure [ 2 ] shows the isosurface of the density function p(x,t ) = \if)(x,t)\ 2 = 0.01 at dif¬ 
ferent times for case 1, while Figure [3] depicts the isosurface evolution for case 2. From 
Figs. |2]j3j we could have the following conclusions: (i). The density of the condensate 
will rotate along with the rotation of the dipole axis. (ii). For Case 1 where the trapping 
potential is isotropic, the shape of the density profile seems to be unchanged during the 
dynamics, and it seems to keep the same symmetric structure with respect to the dipole 
orientation. However, this phenomena does not occur in Case 2 where the trapping po¬ 
tential is anisotropic. 



Figure 2: Isosurface plots of the density function p(x,t ) = \ip(x,t)\ 2 = 0.01 and the dipole 
axis n(f) (red arrow) at different times for case 1 in the example 


4.3 


Example 4.4. Collapse and explosion of a dipolar EEC with damping effect in 3D. 


In this case, the trapping potential V (x) and the constants A and /3 are set to be time 
dependent and are chosen according to the parameters used in the physical experiment 
132 331 (in dimensionless form) as follows: 


V(x,t) = 


(7? * 2 + 7y 3/ 2 + 7- z 2 ) / 2, 

0 , 


t E [0,4 + fhold]/ 

otherwise. 


(4.18) 
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Figure 3: Isosurface plots of the density function p(x,f) = \ ip(x,t)\ 2 = 0.01 and the dipole 
axis n(f) (red arrow) at different times for case 2 in the example 


4.3 


with 7 * = 1.65, 7 y = 1, y z = 1.325, 


where 


A(0 = 


82.864, 

0 , 


t E [0,5.6 +hold]/ 

otherwise. 


0(f) =761.102 < 


1+ 28 

1 1 75£-280' 

tE [0,3.2], 

0.3, 

tE [3.2,3.6], 

1 1 

1 fo(f— 3.6 )' 

tE [3.6,4.8], 

38.8066, 

t E [4.8,5.6 + hoid]/ 

0 , 

otherwise. 


1 f -125t-25e~ 5t +2l5, tE [0,0.4], 

133 | 25 (1 — e“ 2 ) e~ 625t+2 - 5 +140, fe [0.4,1.2]. 


(4.19) 


(4.20) 


(4.21) 


Here, hold is the hold time for the collapse, which is chosen as hold = 0.2. Moreover, we 
let n = (0,0,1) T , n = l and chose the damping term as f(p) = 5p 2 with 5 = 3.512, i.e., we 
chose the cubic nonlinearity and study the case of three-body loss of the particles. The 
initial data in (4.3' is chosen as i/^o( x ) =</ , gs( x ), where tp gs ( x ) is the ground state of the 
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GPE ( |4.1| with f(p ) =0, n = m = n(0), /3 = /3(0) A = A(0) and V (x) = V (x,0), which is 
computed numerically via the numerical method presented in the previous section. The 
computational domain and mesh size are chosen as T>= [—24,24] 3 and h x = h y = /z z = fg, 
respectively. Figures [4] and [5] show the contour plot of the column density 

f R x 

p x c (y r z r t)= \ip(x,t)\ 2 dx, 

J L x 

and the evolution of the total mass, respectively. 



Figure 4: Contour plots of the column density p x c (y,z,t) at different times for the example 
441 


0.9 


0.81 


O 


mass 


2.5 


7.5 


Figure 5: Evolution of the mass for the example 4.4 


From Figs. 4 and 5, we can conclude that: (i). The total mass is lost during the dynam¬ 
ics, especially during a very short period near t = 5 (cf. Fig. 5). (ii). Although the BEC is 
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released from the trap (i.e., the trapping potential is turned off) at time t = 4.2, the atoms 
in the BEC still move inward in the x-y plane, (iii). The density is first enlongated along 
the dipole orientation, then the collapse of the BEC happens very quickly, and "clover" 
pattern of the density profile is created, (iv). The "Leafs" are then ejected outward. All 
these results agree with those in the experiments [32.331. 


5 Conclusion 

We proposed efficient and accurate numerical methods for computing the ground state 
and dynamics of the dipolar Bose-Einstein condensates by integrating a newly devel¬ 
oped dipole-dipole interaction (DDI) solver via the non-uniform fast Fourier transform 
(NUFFT) algorithm [ j28| with existing numerical methods. The NUFFT based DDI solver 
removes naturally the singularity of the DDI at the origin by adopting the spherical/polar 
coordinates in the Fourier space, thus achieves spectral accuracy and simultaneously 
maintains high efficiency by appropriately combining the advantages of the NUFFT and 
FFT. Efficient and accurate numerical methods were then presented to compute the ground 
state and dynamics of the dipolar BEC with a DDI by integrating the normalized gradient 
flow with the backward Euler Fourier pseudospectral discretization and time-splitting 
Fourier pseudospectral method, respectively, together with NUFFT based DDI solver. 
Extensive numerical comparisons with existing methods were carried out to compute 
the DDI, ground states and dynamics of the dipolar BEC. Numerical results showed that 
our new methods outperformed other existing methods in terms of both accuracy and 
efficiency, especially when the computational domain is chosen smaller and/or the solu¬ 
tion is anisotropic. 
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